/* ================================================================================================

matrix.hpp

author: Antoine Allard
affil.: Universite Laval, Quebec City, Quebec
www:    dynamica.phy.ulaval.ca


last modif. August 27, 2013

=================================================================================================*/

#ifndef MATRIX_HPP_INCLUDED
#define MATRIX_HPP_INCLUDED

// standard library
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

using namespace std;


class matrix
{
  protected:
    // dimensions of the matrix
    int dim0, dim1;
    vector<double> dims;
    //
    string mat_filename;
    // stream accessing the file
    ifstream mat_file;
    // string containing a line from the file
    string full_line, value;
    // stream to browse a string
    stringstream one_line;
  protected:
    // function computing the dimensions of the matrix
    void set_dimensions();
    // function filling the elements of the matrix from file
    void fill_matrix(vector< vector< double > >&);
    // function normalizing the elements of the matrix
    void norm_matrix(vector< vector< double > >&);
    // function normalizing each row of the matrix
    void norm_row_matrix(vector< vector< double > >&);
    // adds its transpose to the matrix
    void add_own_transpose(vector< vector< double > >&);
  public:
    // 
    vector<double> get_dimensions();
};

vector<double> matrix::get_dimensions()
{
  vector<double> dims(2);
  dims[0] = dim0;
  dims[1] = dim1;
  return dims;
}

void matrix::set_dimensions()
{
  // initializes the bounds
  dim0 = 0;
  dim1 = 0;

  // stream accessing the file
  ifstream mat_file(mat_filename.c_str());

  // string containing a line from the file
  string full_line, value;

  // stream to browse a string as a stream
  stringstream one_line;

  // reads the file
  if(mat_file.is_open())
  {
    // extracts the first line to compute "j_max"
    getline(mat_file,full_line);
    mat_file >> std::ws;
    ++dim0;

    // computes "j_max"
    one_line.str(full_line);
    while(!one_line.eof())
    {
      one_line >> value >> ws;
      ++dim1;
    }

    // goes through the remaining of the file to compute "i_max"
    while(!mat_file.eof())
    {
      getline(mat_file,full_line);
      mat_file >> ws;
      ++dim0;
    }
    mat_file.close();
  }
  else
  {
    std::cout << "Failed to open file! (" << mat_filename << ")" << std::endl;
    abort();
  }

}

void matrix::fill_matrix(vector< vector< double > >& mat)
{
  //
  mat.resize(dim0);
  //
  mat_file.open(mat_filename.c_str());
  // reads the file
  if(mat_file.is_open())
  {
    //
    for(int i(0); i<dim0; ++i)
    {
      mat[i].resize(dim1);
      getline(mat_file,full_line);
      one_line.str(full_line);
      for(int j(0); j<dim1; ++j)
      {
        one_line >> value >> ws;
        mat[i][j] = atof(value.c_str());
      }
      one_line.clear();
      mat_file >> ws;
    }
    mat_file.close();
  }
  else
  {
    std::cout << "Failed to open file! (" << mat_filename << ")" << std::endl;
    abort();
  }
}

void matrix::norm_matrix(vector< vector< double > >& mat)
{
  double normalisation_value = 0.0;
  for(int i(0); i<dim0; ++i)
    for(int j(0); j<dim1; ++j)
      normalisation_value += mat[i][j];
  for(int i(0); i<dim0; ++i)
    for(int j(0); j<dim1; ++j)
      mat[i][j] /= normalisation_value;
}

void matrix::norm_row_matrix(vector< vector< double > >& mat)
{
  double normalisation_value;
  for(int i(0); i<dim0; ++i)
  {
    normalisation_value = 0.0;
    for(int j(0); j<dim1; ++j)
      normalisation_value += mat[i][j];
    if(normalisation_value > 1e-12)
      for(int j(0); j<dim1; ++j)
        mat[i][j] /= normalisation_value;
  }
}

void matrix::add_own_transpose(vector< vector< double > >& mat)
{
  if(dim0 == dim1)
  {
    for(int i(0); i<dim0; ++i)
    {
      mat[i][i] *= 2.0;
      for(int j(i+1); j<dim1; ++j)
      {
        mat[i][j] += mat[j][i];
        mat[j][i]  = mat[i][j];
      }
    }
  }
}


#endif // MATRIX_HPP_INCLUDED